Setup

# Install packages
# library(remotes)
# remotes::install_github('ropensci/osmdata')

# Data manipulation
# library(tibble)
library(dplyr)
# library(reshape)

# Data visualisation
library(ggplot2)

# Working with spatial data
# library(sf)
# library(leaflet)
library(raster)
library(rgdal)
# library(proj4)
# library(scales)

# Working with OSM data
# library(osmdata)
# library(osmextract)

Intro to raster data (30 + 20 minutes)

View Raster File Attirbutes

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
## rows        1367 
## columns     1697 
## bands       1 
## lower left origin.x        731453 
## lower left origin.y        4712471 
## res.x       1 
## res.y       1 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       1697
## apparent band statistics:
##     Bmin   Bmax    Bmean      Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area
HARV_dsmCrop_info <- capture.output(
  GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
)

Open a Raster in R

We are going to work with a Digital Surface Model (DSM) which is in the GeoTIFF format.

DSM_HARV <- 
  raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")

DSM_HARV
## class      : RasterLayer 
## dimensions : 1367, 1697, 2319799  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_dsmCrop.tif 
## names      : HARV_dsmCrop 
## values     : 305.07, 416.07  (min, max)
summary(DSM_HARV)
##         HARV_dsmCrop
## Min.        305.5500
## 1st Qu.     345.6500
## Median      359.6450
## 3rd Qu.     374.2825
## Max.        413.9000
## NA's          0.0000
summary(DSM_HARV, maxsamp = ncell(DSM_HARV))
##         HARV_dsmCrop
## Min.          305.07
## 1st Qu.       345.59
## Median        359.67
## 3rd Qu.       374.28
## Max.          416.07
## NA's            0.00
DSM_HARV_df <- as.data.frame(DSM_HARV, xy = TRUE)
str(DSM_HARV_df)
## 'data.frame':    2319799 obs. of  3 variables:
##  $ x           : num  731454 731454 731456 731456 731458 ...
##  $ y           : num  4713838 4713838 4713838 4713838 4713838 ...
##  $ HARV_dsmCrop: num  409 408 407 407 409 ...
ggplot() +
    geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = HARV_dsmCrop)) +
    scale_fill_viridis_c() +
    coord_quickmap()

plot(DSM_HARV)

View Raster Coordinate Reference System (CRS) in R

But what units are these?

# CRS have been explained on Day 1
crs(DSM_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 18N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-75,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
##         BBOX[0,-78,84,-72]],
##     ID["EPSG",32618]]

Calculate the Min and Max value

minValue(DSM_HARV)
## [1] 305.07
maxValue(DSM_HARV)
## [1] 416.07

If Min and Max values haven’t been calculated, you can set them with the raster::setMinMax() function.

# Maybe skip this
DSM_HARV <- raster::setMinMax(DSM_HARV)

Raster bands

To see how many bands a raster dataset has, use the raster::nlayers() function.

nlayers(DSM_HARV)
## [1] 1

This dataset has only 1 band. We will discuss multi-band raster data in a later episode.

Dealing with missing data

In raster data, pixels with no value are represented with NoDataValue. That is usually the case with edge values where data is missing from the rectangular region of the raster. Use raster::GDALinfo() to check the metadata for missing data values.

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
## rows        1367 
## columns     1697 
## bands       1 
## lower left origin.x        731453 
## lower left origin.y        4712471 
## res.x       1 
## res.y       1 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       1697
## apparent band statistics:
##     Bmin   Bmax    Bmean      Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area

Bad data values

Bad values are usually miscalculations that result in values out of a predefined range.

# This code chunk is hidden in the lesson material

Creating a histogram of raster values

A histogram can be used to inspect the distribution of raster values visually. It can show if there are values above the max or below the min of the expected range. We can create a histogram with the ggplot2 function geom_histogram().

ggplot() +
    geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop))

Adjust the level of desired detail by setting the number of bins.

ggplot() +
    geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop), bins = 40)

Challenge

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
## rows        1367 
## columns     1697 
## bands       1 
## lower left origin.x        731453 
## lower left origin.y        4712471 
## res.x       1 
## res.y       1 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       1697
## apparent band statistics:
##         Bmin      Bmax     Bmean       Bsd
## 1 -0.7136298 0.9999997 0.3125525 0.4812939
## Metadata:
## AREA_OR_POINT=Area

Plotting raster data (40 + 30 minutes - 10 minutes about plot formatting - 15 minutes details = 45 minutes)

DSM_HARV_df <- DSM_HARV_df %>%
                mutate(fct_elevation = cut(HARV_dsmCrop, breaks = 3))

ggplot() +
    geom_bar(data = DSM_HARV_df, aes(fct_elevation))

To see the cutoff values:

unique(DSM_HARV_df$fct_elevation)
## [1] (379,416] (342,379] (305,342]
## Levels: (305,342] (342,379] (379,416]

To show count number of pixels in each group:

DSM_HARV_df %>%
        group_by(fct_elevation) %>%
        count()
## # A tibble: 3 × 2
## # Groups:   fct_elevation [3]
##   fct_elevation       n
##   <fct>           <int>
## 1 (305,342]      418891
## 2 (342,379]     1530073
## 3 (379,416]      370835

To customize cutoff values:

custom_bins <- c(300, 350, 400, 450)

DSM_HARV_df <- DSM_HARV_df %>%
  mutate(fct_elevation_2 = cut(HARV_dsmCrop, breaks = custom_bins))

unique(DSM_HARV_df$fct_elevation_2)
## [1] (400,450] (350,400] (300,350]
## Levels: (300,350] (350,400] (400,450]
ggplot() +
  geom_bar(data = DSM_HARV_df, aes(fct_elevation_2))

DSM_HARV_df %>%
  group_by(fct_elevation_2) %>%
  count()
## # A tibble: 3 × 2
## # Groups:   fct_elevation_2 [3]
##   fct_elevation_2       n
##   <fct>             <int>
## 1 (300,350]        741815
## 2 (350,400]       1567316
## 3 (400,450]         10668
ggplot() +
  geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = fct_elevation_2)) + 
  coord_quickmap()

Customising colors:

terrain.colors(3)
## [1] "#00A600" "#ECB176" "#F2F2F2"
ggplot() +
 geom_raster(data = DSM_HARV_df , aes(x = x, y = y,
                                      fill = fct_elevation_2)) + 
    scale_fill_manual(values = terrain.colors(3)) + 
    coord_quickmap()

my_col <- terrain.colors(3)

ggplot() +
 geom_raster(data = DSM_HARV_df , aes(x = x, y = y,
                                      fill = fct_elevation_2)) + 
    scale_fill_manual(values = my_col, name = "Elevation") + 
    coord_quickmap()

Challenge

DSM_HARV_df <- DSM_HARV_df %>%
  mutate(fct_elevation_6 = cut(HARV_dsmCrop, breaks = 6))

unique(DSM_HARV_df$fct_elevation_6)
## [1] (398,416] (379,398] (361,379] (342,361] (324,342] (305,324]
## Levels: (305,324] (324,342] (342,361] (361,379] (379,398] (398,416]
my_col <- terrain.colors(6)

ggplot() +
  geom_raster(data = DSM_HARV_df, aes(x = x, y = y,
                                       fill = fct_elevation_6)) +
  scale_fill_manual(values = my_col, name = "Elevation") +
  coord_quickmap() +
  xlab("X") +
  ylab("Y") +
  labs(title = "Elevation Classes of the Harvard Forest Digital Surface Model (DSM)")

Layering rasters

DSM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")

DSM_hill_HARV
## class      : RasterLayer 
## dimensions : 1367, 1697, 2319799  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_DSMhill.tif 
## names      : HARV_DSMhill 
## values     : -0.7136298, 0.9999997  (min, max)
DSM_hill_HARV_df <- as.data.frame(DSM_hill_HARV, xy = TRUE)

str(DSM_hill_HARV_df)
## 'data.frame':    2319799 obs. of  3 variables:
##  $ x           : num  731454 731454 731456 731456 731458 ...
##  $ y           : num  4713838 4713838 4713838 4713838 4713838 ...
##  $ HARV_DSMhill: num  NA NA NA NA NA NA NA NA NA NA ...
ggplot() +
  geom_raster(data = DSM_hill_HARV_df,
              aes(x = x, y = y, alpha = HARV_DSMhill)) + 
  scale_alpha(range =  c(0.15, 0.65), guide = "none") + 
  coord_quickmap()

ggplot() +
  geom_raster(data = DSM_HARV_df , 
              aes(x = x, y = y, 
                  fill = HARV_dsmCrop)) + 
  geom_raster(data = DSM_hill_HARV_df, 
              aes(x = x, y = y, 
                  alpha = HARV_DSMhill)) +  
  scale_fill_viridis_c() +  
  scale_alpha(range = c(0.15, 0.65), guide = "none") +  
  ggtitle("Elevation with hillshade") +
  coord_quickmap()

Challenge

# CREATE DSM MAPS

# import DSM data
DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
# convert to a df for plotting
DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)

# import DSM hillshade
DSM_hill_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmHill.tif")
# convert to a df for plotting
DSM_hill_SJER_df <- as.data.frame(DSM_hill_SJER, xy = TRUE)

# Build Plot
ggplot() +
    geom_raster(data = DSM_SJER_df , 
                aes(x = x, y = y, 
                     fill = SJER_dsmCrop,
                     alpha = 0.8)
                ) + 
    geom_raster(data = DSM_hill_SJER_df, 
                aes(x = x, y = y, 
                  alpha = SJER_dsmHill)
                ) +
    scale_fill_viridis_c() +
    guides(fill = guide_colorbar()) +
    scale_alpha(range = c(0.4, 0.7), guide = "none") +
    # remove grey background and grid lines
    theme_bw() + 
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank()) +
    xlab("UTM Easting Coordinate (m)") +
    ylab("UTM Northing Coordinate (m)") +
    ggtitle("DSM with Hillshade") +
    coord_quickmap()

# CREATE DTM MAP
# import DTM
DTM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmCrop.tif")
DTM_SJER_df <- as.data.frame(DTM_SJER, xy = TRUE)

# DTM Hillshade
DTM_hill_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmHill.tif")
DTM_hill_SJER_df <- as.data.frame(DTM_hill_SJER, xy = TRUE)

ggplot() +
    geom_raster(data = DTM_SJER_df ,
                aes(x = x, y = y,
                     fill = SJER_dtmCrop,
                     alpha = 2.0)
                ) +
    geom_raster(data = DTM_hill_SJER_df,
                aes(x = x, y = y,
                  alpha = SJER_dtmHill)
                ) +
    scale_fill_viridis_c() +
    guides(fill = guide_colorbar()) +
    scale_alpha(range = c(0.4, 0.7), guide = "none") +
    theme_bw() +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank()) +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank()) +
    ggtitle("DTM with Hillshade") +
    coord_quickmap()

Reproject raster data (40 + 20 minutes - 10 minutes instruction details - 20 minutes challenges = 30 minutes)

What happens when maps don’t line up? That is usually a sign that layers are in different coordinate reference systems (CRS).

DTM_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")

DTM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")
DTM_HARV_df <- as.data.frame(DTM_HARV, xy = TRUE)

DTM_hill_HARV_df <- as.data.frame(DTM_hill_HARV, xy = TRUE)
ggplot() +
     geom_raster(data = DTM_HARV_df , 
                 aes(x = x, y = y, 
                  fill = HARV_dtmCrop)) + 
     geom_raster(data = DTM_hill_HARV_df, 
                 aes(x = x, y = y, 
                   alpha = HARV_DTMhill_WGS84)) +
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
     coord_quickmap()

ggplot() +
geom_raster(data = DTM_HARV_df,
    aes(x = x, y = y,
    fill = HARV_dtmCrop)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
coord_quickmap()

ggplot() +
geom_raster(data = DTM_hill_HARV_df,
    aes(x = x, y = y,
    alpha = HARV_DTMhill_WGS84)) + 
    coord_quickmap()

### Challenge

crs(DTM_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 18N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-75,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
##         BBOX[0,-78,84,-72]],
##     ID["EPSG",32618]]
crs(DTM_hill_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##     REMARK["Axis order reversed compared to EPSG:4326"]]

Reproject rasters

DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
                                       crs = crs(DTM_HARV))
crs(DTM_hill_UTMZ18N_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 18N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-75,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
##         BBOX[0,-78,84,-72]],
##     ID["EPSG",32618]]
crs(DTM_hill_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##     REMARK["Axis order reversed compared to EPSG:4326"]]
extent(DTM_hill_UTMZ18N_HARV)
## class      : Extent 
## xmin       : 731397.3 
## xmax       : 733205.3 
## ymin       : 4712403 
## ymax       : 4713907
extent(DTM_hill_HARV)
## class      : Extent 
## xmin       : -72.18192 
## xmax       : -72.16061 
## ymin       : 42.52941 
## ymax       : 42.54234

Challenge

Dealing with raster resolution

res(DTM_hill_UTMZ18N_HARV)
## [1] 1.000 0.998
res(DTM_HARV)
## [1] 1 1
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
                                         crs = crs(DTM_HARV),
                                         res = res(DTM_HARV))
res(DTM_hill_UTMZ18N_HARV)
## [1] 1 1
res(DTM_HARV)
## [1] 1 1
DTM_hill_HARV_2_df <- as.data.frame(DTM_hill_UTMZ18N_HARV, xy = TRUE)
ggplot() +
     geom_raster(data = DTM_HARV_df , 
                 aes(x = x, y = y, 
                  fill = HARV_dtmCrop)) + 
     geom_raster(data = DTM_hill_HARV_2_df, 
                 aes(x = x, y = y, 
                   alpha = HARV_DTMhill_WGS84)) +
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
     coord_quickmap()

Challenge

# import DSM
DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
# import DSM hillshade
DSM_hill_SJER_WGS <-
raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")

# reproject raster
DTM_hill_UTMZ18N_SJER <- projectRaster(DSM_hill_SJER_WGS,
                                  crs = crs(DSM_SJER),
                                  res = 1)

# convert to data.frames
DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)

DSM_hill_SJER_df <- as.data.frame(DTM_hill_UTMZ18N_SJER, xy = TRUE)

ggplot() +
     geom_raster(data = DSM_hill_SJER_df, 
                 aes(x = x, y = y, 
                   alpha = SJER_DSMhill_WGS84)
                 ) +
     geom_raster(data = DSM_SJER_df, 
             aes(x = x, y = y, 
                  fill = SJER_dsmCrop,
                  alpha=0.8)
             ) + 
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
     coord_quickmap()

Raster calculations (40 + 20 minutes - 15 minutes raster math - 15 minutes challenges = 30 minutes)

Raster calculations in R

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
## rows        1367 
## columns     1697 
## bands       1 
## lower left origin.x        731453 
## lower left origin.y        4712471 
## res.x       1 
## res.y       1 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       1697
## apparent band statistics:
##     Bmin   Bmax    Bmean      Bsd
## 1 304.56 389.82 344.8979 15.86147
## Metadata:
## AREA_OR_POINT=Area
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
## rows        1367 
## columns     1697 
## bands       1 
## lower left origin.x        731453 
## lower left origin.y        4712471 
## res.x       1 
## res.y       1 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       1697
## apparent band statistics:
##     Bmin   Bmax    Bmean      Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area
 ggplot() +
      geom_raster(data = DTM_HARV_df , 
              aes(x = x, y = y, fill = HARV_dtmCrop)) +
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
     coord_quickmap()

 ggplot() +
      geom_raster(data = DSM_HARV_df , 
              aes(x = x, y = y, fill = HARV_dsmCrop)) +
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
     coord_quickmap()

Raster math and Canopy Height Models (DROP THIS PART)

CHM_HARV <- DSM_HARV - DTM_HARV

CHM_HARV_df <- as.data.frame(CHM_HARV, xy = TRUE)
 ggplot() +
   geom_raster(data = CHM_HARV_df , 
               aes(x = x, y = y, fill = layer)) + 
   scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
   coord_quickmap()

ggplot(CHM_HARV_df) +
    geom_histogram(aes(layer))

Challenge

custom_bins <- c(0, 10, 20, 30, 40)
CHM_HARV_df <- CHM_HARV_df %>%
                  mutate(canopy_discrete = cut(layer, breaks = custom_bins))

ggplot() +
  geom_raster(data = CHM_HARV_df , aes(x = x, y = y,
                                       fill = canopy_discrete)) + 
     scale_fill_manual(values = terrain.colors(4)) + 
     coord_quickmap()

Efficient raster calculations: overlay function

CHM_ov_HARV <- overlay(DSM_HARV,
                       DTM_HARV,
                       fun = function(r1, r2) { return( r1 - r2) })
CHM_ov_HARV_df <- as.data.frame(CHM_ov_HARV, xy = TRUE)
 ggplot() +
   geom_raster(data = CHM_ov_HARV_df, 
               aes(x = x, y = y, fill = layer)) + 
   scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
   coord_quickmap()

Export a GeoTIFF

writeRaster(CHM_ov_HARV, "CHM_HARV.tiff",
            format="GTiff",
            overwrite=TRUE,
            NAflag=-9999)

Challenge

CHM_ov_SJER <- overlay(DSM_SJER,
                       DTM_SJER,
                       fun = function(r1, r2){ return(r1 - r2) })
CHM_ov_SJER_df <- as.data.frame(CHM_ov_SJER, xy = TRUE)
ggplot(CHM_ov_SJER_df) +
    geom_histogram(aes(layer))

ggplot() +
  geom_raster(data = CHM_ov_SJER_df,
              aes(x = x, y = y,
                  fill = layer)) +
  scale_fill_gradientn(name = "Canopy Height",
                       colors = terrain.colors(10)) +
  coord_quickmap()

writeRaster(CHM_ov_SJER, "chm_ov_SJER.tiff",
            format = "GTiff",
            overwrite = TRUE,
            NAflag = -9999)
ggplot(CHM_HARV_df) +
    geom_histogram(aes(layer))

ggplot(CHM_ov_SJER_df) +
    geom_histogram(aes(layer))

Work with multi-band rasters (40 + 20 minutes - 15 minutes details - 15 minutes challenges = 30 minutes)

##Getting Started with Multi-Band Data in R

The raster() function only reads in the first band, in this case the red band of an RGB raster

RGB_band1_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
RGB_band1_HARV_df  <- as.data.frame(RGB_band1_HARV, xy = TRUE)
ggplot() +
  geom_raster(data = RGB_band1_HARV_df,
              aes(x = x, y = y, alpha = HARV_RGB_Ortho)) + 
  coord_quickmap()

RGB_band1_HARV
## class      : RasterLayer 
## band       : 1  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho 
## values     : 0, 255  (min, max)
nbands(RGB_band1_HARV)
## [1] 3

To import the green band:

RGB_band2_HARV <-  raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif", band = 2)
RGB_band2_HARV_df <- as.data.frame(RGB_band2_HARV, xy = TRUE)
ggplot() +
  geom_raster(data = RGB_band2_HARV_df,
              aes(x = x, y = y, alpha = HARV_RGB_Ortho)) + 
  coord_equal()

Raster stacks

The stack() function brings in all bands

RGB_stack_HARV <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
RGB_stack_HARV
## class      : RasterStack 
## dimensions : 2317, 3073, 7120141, 3  (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## names      : HARV_RGB_Ortho.1, HARV_RGB_Ortho.2, HARV_RGB_Ortho.3 
## min values :                0,                0,                0 
## max values :              255,              255,              255
RGB_stack_HARV@layers
## [[1]]
## class      : RasterLayer 
## band       : 1  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.1 
## values     : 0, 255  (min, max)
## 
## 
## [[2]]
## class      : RasterLayer 
## band       : 2  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.2 
## values     : 0, 255  (min, max)
## 
## 
## [[3]]
## class      : RasterLayer 
## band       : 3  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.3 
## values     : 0, 255  (min, max)
RGB_stack_HARV[[2]]
## class      : RasterLayer 
## band       : 2  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : HARV_RGB_Ortho.tif 
## names      : HARV_RGB_Ortho.2 
## values     : 0, 255  (min, max)
RGB_stack_HARV_df  <- as.data.frame(RGB_stack_HARV, xy = TRUE)
str(RGB_stack_HARV_df)
## 'data.frame':    7120141 obs. of  5 variables:
##  $ x               : num  731999 731999 731999 731999 732000 ...
##  $ y               : num  4713535 4713535 4713535 4713535 4713535 ...
##  $ HARV_RGB_Ortho.1: num  0 2 6 0 16 0 0 6 1 5 ...
##  $ HARV_RGB_Ortho.2: num  1 0 9 0 5 0 4 2 1 0 ...
##  $ HARV_RGB_Ortho.3: num  0 10 1 0 17 0 4 0 0 7 ...
ggplot() +
  geom_histogram(data = RGB_stack_HARV_df, aes(HARV_RGB_Ortho.1))

ggplot() +
  geom_raster(data = RGB_stack_HARV_df,
              aes(x = x, y = y, alpha = HARV_RGB_Ortho.2)) + 
  coord_quickmap()

Create a three-band image

plotRGB(RGB_stack_HARV,
        r = 1, g = 2, b = 3)

Which is the same with terra

terra::plotRGB(RGB_stack_HARV,
        r = 1, g = 2, b = 3)

plotRGB(RGB_stack_HARV,
        r = 1, g = 2, b = 3,
        scale = 800,
        stretch = "lin")

plotRGB(RGB_stack_HARV,
        r = 1, g = 2, b = 3,
        scale = 800,
        stretch = "hist")

Challenge

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")
## rows        2317 
## columns     3073 
## bands       3 
## lower left origin.x        731998.5 
## lower left origin.y        4712956 
## res.x       0.25 
## res.y       0.25 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE       -9999          1       3073
## 2 Float64           TRUE       -9999          1       3073
## 3 Float64           TRUE       -9999          1       3073
## apparent band statistics:
##   Bmin Bmax     Bmean      Bsd
## 1    0  255 107.83651 30.01918
## 2    0  255 130.09595 32.00168
## 3    0  255  95.75979 16.57704
## Metadata:
## AREA_OR_POINT=Area
HARV_NA <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")
plotRGB(HARV_NA,
        r = 1, g = 2, b = 3)

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
## rows        2317 
## columns     3073 
## bands       3 
## lower left origin.x        731998.5 
## lower left origin.y        4712956 
## res.x       0.25 
## res.y       0.25 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE   -1.7e+308          1       3073
## 2 Float64           TRUE   -1.7e+308          1       3073
## 3 Float64           TRUE   -1.7e+308          1       3073
## apparent band statistics:
##   Bmin Bmax Bmean Bsd
## 1    0  255   NaN NaN
## 2    0  255   NaN NaN
## 3    0  255   NaN NaN
## Metadata:
## AREA_OR_POINT=Area

RasterStack vs. RasterBrick

object.size(RGB_stack_HARV)
## 51296 bytes
RGB_brick_HARV <- brick(RGB_stack_HARV)

object.size(RGB_brick_HARV)
## 170898912 bytes
plotRGB(RGB_brick_HARV)

methods(class=class(RGB_stack_HARV))
##   [1] !                     !=                    [                    
##   [4] [[                    [[<-                  [<-                  
##   [7] %in%                  ==                    $                    
##  [10] $<-                   addLayer              adjacent             
##  [13] aggregate             all.equal             animate              
##  [16] approxNA              area                  Arith                
##  [19] as.array              as.character          as.data.frame        
##  [22] as.integer            as.list               as.logical           
##  [25] as.matrix             as.vector             atan2                
##  [28] bbox                  blockSize             boxplot              
##  [31] brick                 calc                  cellFromRowCol       
##  [34] cellFromRowColCombine cellFromXY            cellStats            
##  [37] clamp                 click                 coerce               
##  [40] colFromCell           colFromX              colSums              
##  [43] Compare               coordinates           corLocal             
##  [46] couldBeLonLat         cover                 crop                 
##  [49] crosstab              crs<-                 cut                  
##  [52] cv                    density               dim                  
##  [55] dim<-                 disaggregate          dropLayer            
##  [58] extend                extent                extract              
##  [61] flip                  freq                  getValues            
##  [64] getValuesBlock        getValuesFocal        hasValues            
##  [67] head                  hist                  image                
##  [70] init                  inMemory              interpolate          
##  [73] intersect             is.factor             is.finite            
##  [76] is.infinite           is.na                 is.nan               
##  [79] isLonLat              KML                   labels               
##  [82] length                levels                levels<-             
##  [85] log                   Logic                 mask                 
##  [88] match                 Math                  Math2                
##  [91] maxValue              mean                  merge                
##  [94] metadata              minValue              modal                
##  [97] mosaic                names                 names<-              
## [100] ncell                 ncol                  ncol<-               
## [103] nlayers               nrow                  nrow<-               
## [106] origin                origin<-              overlay              
## [109] pairs                 persp                 plot                 
## [112] plotRGB               predict               print                
## [115] proj4string           proj4string<-         quantile             
## [118] raster                rasterize             ratify               
## [121] readAll               readStart             readStop             
## [124] reclassify            rectify               res                  
## [127] res<-                 resample              rotate               
## [130] rowColFromCell        rowFromCell           rowFromY             
## [133] rowSums               sampleRandom          sampleRegular        
## [136] scale                 select                setMinMax            
## [139] setValues             shift                 show                 
## [142] spplot                stack                 stackSelect          
## [145] stretch               subs                  subset               
## [148] Summary               summary               t                    
## [151] tail                  text                  trim                 
## [154] unique                unstack               values               
## [157] values<-              weighted.mean         which.max            
## [160] which.min             whiches.max           whiches.min          
## [163] wkt                   writeRaster           xFromCell            
## [166] xFromCol              xmax                  xmax<-               
## [169] xmin                  xmin<-                xres                 
## [172] xyFromCell            yFromCell             yFromRow             
## [175] ymax                  ymax<-                ymin                 
## [178] ymin<-                yres                  zonal                
## [181] zoom                 
## see '?methods' for accessing help and source code
methods(class=class(RGB_stack_HARV[1]))
##  [1] [             [<-           anyDuplicated as_tibble     as.data.frame
##  [6] as.raster     bind          boxplot       brick         cellFromXY   
## [11] coerce        coordinates   determinant   distance      duplicated   
## [16] edit          extent        extract       head          initialize   
## [21] isSymmetric   Math          Math2         Ops           raster       
## [26] rasterize     relist        subset        summary       surfaceArea  
## [31] tail          trim          unique        values<-      weighted.mean
## [36] writeValues  
## see '?methods' for accessing help and source code